%Script to produce publication-quality figures for Kiefer et al. (2025)
%manuscript.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PRELIMINARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
b=pwd;b=strsplit(b,'/');b=b{2};
if(strcmp(b,'home'))%102 server (CPU and memory-intensive processing and/or figure generation)
    addpath(genpath('/home/mmfire/matlab_functions/FUNCTIONS'))
    addpath(genpath('/home/mmfire/matlab_functions/SHAPEFILES'))
    arps_dir='/d102c/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL_output/ARPS-CANOPY/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/';
    devs_dir='/d102c/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL_output/DEVS-FIRE/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/';
    grid_dir='/home/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-DEVS_IDEALIZED_TESTS/PHASE2/GRIDINFO_ANALYSIS/';%GRID INFO IS IDENTICAL FOR ARPS and ARPS-CANOPY SIMS.
    plot_dir='/home/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS_CODEDEVEL/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL/diff_plots/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/COUPLING_ANALYSIS/SURF_tightview/';
elseif(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server (CPU and memory-intensive processing and/or figure generation)
    addpath(genpath('/dska/mtkiefer/matlab_functions/FUNCTIONS'))
    addpath(genpath('/dska/mtkiefer/matlab_functions/SHAPEFILES'))
    arps_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS_output/SERDP_2019_ACDF/ARPS-CANOPY/RUNS/TEST4/';
    devs_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS_output/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    devsconfig_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    grid_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/GRID_INFO/';
    plot_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/PLOTS/TEST4/';
    mat_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/MATFILES/TEST4/';
elseif(strcmp(b,'Volumes'))%my laptop (simpler processing and/or figure generation)
    addpath(genpath('/Users/mmfire/Documents/MATLAB/matlab_functions/FUNCTIONS'))
    addpath(genpath('/Users/mmfire/Documents/MATLAB/matlab_functions/SHAPEFILES'))
    arps_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/RUNS/ARPS-CANOPY/TEST4/';
    devs_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/RUNS/DEVS-FIRE/TEST4/';
    devsconfig_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    grid_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/GRID_INFO/';
    plot_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/PLOTS/TEST4/';
    mat_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/MATFILES/TEST4/';
end
arps_head='acdf_serdp_0030';
%Check if initial processing has already been performed.
if(~isfile(strcat(mat_dir,'/forpub/acfire_figures_forpub.mat'))) %REQUISITE MAT FILE NOT FOUND
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %1. Read in grid coordinate information.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ARPS.DX=30;%ARPS horizontal grid spacing
    DEVS.DX=3;%DEVS-FIRE horizontal grid spacing
    RESRATX=ARPS.DX/DEVS.DX;RESRATY=RESRATX;%ratio of ARPS to DEVS-FIRE horizontal grid spacing.
    NX=253;NY=253;NZ=83;%ARPS grid dimensions
    NWRITEPTS_X=500;NWRITEPTS_Y=500;%DEVS-FIRE grid dimensions
    %1a. ARPS
    %Note: Coordinates are valid at the center of each grid cell, in
    %conformance with contourm standard.
    load('/dska/mtkiefer/SERDP/CASES/SLEF_2019/PROD_ANALYSIS/FILES_MAT/plot_domainmaps_forpub.mat')
    ARPS.LATS=D5.XLATF;ARPS.LATS(253,253)=NaN;%i=NX,j=NY lopped off by ARPS when dumping out to NetCDF files.
    ARPS.LONS=D5.XLONGF;ARPS.LONS(253,253)=NaN;%i=NX,j=NY lopped off by ARPS when dumping out to NetCDF files.
    clear D1 D2 D3 D4 D5
    ARPS.X=nan(size(ARPS.LATS,1),size(ARPS.LATS,2));ARPS.Y=nan(size(ARPS.LATS,1),size(ARPS.LATS,2));
    for j=1:size(ARPS.LATS,1)
        ARPS.X(j,:)=[-ARPS.DX:ARPS.DX:ARPS.DX*size(ARPS.X,2)-(2*ARPS.DX)];%LL corner of cell
    end
    for i=1:size(ARPS.LATS,2)
        ARPS.Y(:,i)=[-ARPS.DX:ARPS.DX:ARPS.DX*size(ARPS.X,2)-(2*ARPS.DX)]';%LL corner of cell
    end
    %1b. DEVS
    %Note: Coordinates are valid at the lower-left corner of each grid cell,
    %in conformance with pcolorm standard.
    filename = strcat(grid_dir,'gridinfo_devs.txt');
    formatSpec = '%5f%8f%11f%11f%15f%f%[^\n\r]';
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'EmptyValue', NaN,  'ReturnOnError', false);
    fclose(fileID);
    temW = [dataArray{1:end-1}];
    clearvars filename formatSpec fileID dataArray ans;
    DEVS.X=nan(NWRITEPTS_Y+6*RESRATY,NWRITEPTS_X+4*RESRATX);DEVS.Y=nan(NWRITEPTS_Y+6*RESRATY,NWRITEPTS_X+4*RESRATX);%Source: arps.f90
    DEVS.LONS=nan(NWRITEPTS_Y+6*RESRATY,NWRITEPTS_X+4*RESRATX);DEVS.LATS=nan(NWRITEPTS_Y+6*RESRATY,NWRITEPTS_X+4*RESRATX);%Source: arps.f90
    for r=1:size(temW,1)
        i=temW(r,1);j=temW(r,2);
        DEVS.X(j,i)=temW(r,3);
        DEVS.Y(j,i)=temW(r,4);
        DEVS.LONS(j,i)=temW(r,5);
        DEVS.LATS(j,i)=temW(r,6);
        clear i j
    end
    clear r tem
    %ARPS pads the DEVS grid coordinate arrays. Here I remove the extra
    %elements. These index ranges place the ignition line in the *identical*
    %position in the ARPS physical domain as that derived from the surface heat
    %flux variable in the ARPS history files.
    DEVS.X=DEVS.X(2*RESRATY+2:NWRITEPTS_Y+2*RESRATY+1,2*RESRATX+2:NWRITEPTS_X+2*RESRATX+1);
    DEVS.Y=DEVS.Y(2*RESRATY+2:NWRITEPTS_Y+2*RESRATY+1,2*RESRATX+2:NWRITEPTS_X+2*RESRATX+1);
    DEVS.LONS=DEVS.LONS(2*RESRATY+2:NWRITEPTS_Y+2*RESRATY+1,2*RESRATX+2:NWRITEPTS_X+2*RESRATX+1);
    DEVS.LATS=DEVS.LATS(2*RESRATY+2:NWRITEPTS_Y+2*RESRATY+1,2*RESRATX+2:NWRITEPTS_X+2*RESRATX+1);
    %%%PAUSE TO MAKE SURE THAT GRID COORDINATE DEFINITIONS ARE CONSISTENT.  I
    %%%WANT LAT2 AND LON2 TO BE THE COORDINATES AT THE *LOWER LEFT* CORNER OF
    %%%THE GRID CELLS.
    DEVS.LAT2=DEVS.LATS;DEVS.LON2=DEVS.LONS;
    DEVS.LATS=nan(1)*DEVS.LAT2;for j=1:size(DEVS.LATS,1)-1;for i=1:size(DEVS.LATS,2)-1;DEVS.LATS(j,i)=0.25*(DEVS.LAT2(j,i)+DEVS.LAT2(j+1,i)+DEVS.LAT2(j,i+1)+DEVS.LAT2(j+1,i+1));end;end
    DEVS.LONS=nan(1)*DEVS.LON2;for j=1:size(DEVS.LONS,1)-1;for i=1:size(DEVS.LONS,2)-1;DEVS.LONS(j,i)=0.25*(DEVS.LON2(j,i)+DEVS.LON2(j+1,i)+DEVS.LON2(j,i+1)+DEVS.LON2(j+1,i+1));end;end
    %For pcolorm, the lower-left corner coordinates are correct.  For
    %quivermc, the coordinates should be at the center of the grid cell.
    ARPS.LAT2=nan(1)*ARPS.LATS;for j=2:size(ARPS.LAT2,1);for i=2:size(ARPS.LAT2,2);ARPS.LAT2(j,i)=0.25*(ARPS.LATS(j,i)+ARPS.LATS(j-1,i)+ARPS.LATS(j,i-1)+ARPS.LATS(j-1,i-1));end;end
    ARPS.LON2=nan(1)*ARPS.LONS;for j=2:size(ARPS.LON2,1);for i=2:size(ARPS.LON2,2);ARPS.LON2(j,i)=0.25*(ARPS.LONS(j,i)+ARPS.LONS(j-1,i)+ARPS.LONS(j,i-1)+ARPS.LONS(j-1,i-1));end;end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %2. Read in heatmap data from DEVS-FIRE heatmap files every second: DEVS.BPRG
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ts=0;ti=1;te=7200;
    SHFX=nan(NWRITEPTS_Y,NWRITEPTS_X,(te-ts)/ti+1);
    c=2;%skip t=0 (empty heatmap file)
    for t=ts+ti:ti:te%skip t=0 (empty heatmap file)
        disp(num2str(t))
        if(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server
            filename = strcat(devs_dir,num2str(t),'.txt');
            delimiter = '\t';
            startRow = 7;
            formatSpec = '%f%f%f%f%f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            textscan(fileID, '%[^\n\r]', startRow-1, 'WhiteSpace', '', 'ReturnOnError', false, 'EndOfLine', '\r\n');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false);
            fclose(fileID);
            temW = [dataArray{1:end-1}];
            clearvars filename delimiter startRow formatSpec fileID dataArray ans;
        elseif(strcmp(b,'Volumes'))%my laptop
            opts = delimitedTextImportOptions("NumVariables", 6);
            opts.DataLines = [8, Inf];
            opts.Delimiter = "\t";
            opts.VariableNames = ["CellSideSizem", "SpaceWidth", "SpaceHeight", "VarName4", "VarName5", "VarName6"];
            opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];
            opts.ExtraColumnsRule = "ignore";
            opts.EmptyLineRule = "read";
            temW = readtable(strcat(devs_dir,num2str(t),'.txt'), opts);
            temW = table2array(temW);
            clear opts
        end
        for r=1:size(temW,1)
            if(~isnan(temW(r,2)))
                SHFX(temW(r,2),temW(r,1),c)=temW(r,3);
            end
        end
        c=c+1;
    end
    %Compute array containing the time (timesteps since initialization) when
    %the sensible heat flux in each DEVS grid cell first exceeds 500 W/m2
    DEVS.BPRG=nan(size(SHFX,1),size(SHFX,2));%Burn progress
    for j=1:size(SHFX,1)
        for i=1:size(SHFX,2)
            idx = find(SHFX(j,i,:)>500,1,'first');
            if(~isempty(idx)==1)
                DEVS.BPRG(j,i)=idx;
            else
                DEVS.BPRG(j,i)=9e99;%MTK 7/16/24: 9e99 = never-burn cells.
            end
        end
    end
    clear SHFX j i idx c ts ti te temW
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %4. Read in heatmap data from DEVS-FIRE heatmap files every minute: DEVS.HS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ts=0;ti=60;te=7200;%Read in every minute.
    DEVS.HS=nan(NWRITEPTS_Y,NWRITEPTS_X,(te-ts)/ti+1);
    c=2;%skip t=0 (empty heatmap file)
    for t=ts+ti:ti:te%skip t=0 (empty heatmap file)
        if(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server
            filename = strcat(devs_dir,num2str(t),'.txt');
            delimiter = '\t';
            startRow = 7;
            formatSpec = '%f%f%f%f%f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            textscan(fileID, '%[^\n\r]', startRow-1, 'WhiteSpace', '', 'ReturnOnError', false, 'EndOfLine', '\r\n');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false);
            fclose(fileID);
            temW = [dataArray{1:end-1}];
            clearvars filename delimiter startRow formatSpec fileID dataArray ans;
        elseif(strcmp(b,'Volumes'))%my laptop
            opts = delimitedTextImportOptions("NumVariables", 6);
            opts.DataLines = [8, Inf];
            opts.Delimiter = "\t";
            opts.VariableNames = ["CellSideSizem", "SpaceWidth", "SpaceHeight", "VarName4", "VarName5", "VarName6"];
            opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];
            opts.ExtraColumnsRule = "ignore";
            opts.EmptyLineRule = "read";
            temW = readtable(strcat(devs_dir,num2str(t),'.txt'), opts);
            temW = table2array(temW);
            clear opts
        end
        for r=1:size(temW,1)
            if(~isnan(temW(r,2)))
                DEVS.HS(temW(r,2),temW(r,1),c)=temW(r,3);
            end
        end
        c=c+1;
    end
    %%%Mask values less than 500 W/m2
    DEVS.HS(DEVS.HS<500)=NaN;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %5. Read in interpolated winds from weather.txt files output by ARPS-CANOPY
    %(DEVS-FIRE grid cells): DEVS.U6, DEVS.V6
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ts=0;ti=60;te=7200;%Read in every minute.
    DEVS.WSPD=nan(NWRITEPTS_Y,NWRITEPTS_X,(te-ts)/ti+1);
    DEVS.WDIR=nan(NWRITEPTS_Y,NWRITEPTS_X,(te-ts)/ti+1);
    c=2;%skip t=0 (empty heatmap file)
    for t=ts+ti:ti:te%skip t=0 (empty heatmap file)
        disp(sprintf('%d',t))
        if(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server
            opts = delimitedTextImportOptions("NumVariables", 10);
            opts.DataLines = [2, Inf];
            opts.Delimiter = " ";
            opts.VariableNames = ["X_PT", "Y_PT", "LATITUDE", "LONGITUDE", "TEMPA", "TEMPS", "RELH", "WSPD", "WDIR", "Var10"];
            opts.SelectedVariableNames = ["X_PT", "Y_PT", "LATITUDE", "LONGITUDE", "TEMPA", "TEMPS", "RELH", "WSPD", "WDIR"];
            opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "string"];
            opts.ExtraColumnsRule = "ignore";
            opts.EmptyLineRule = "read";
            opts.ConsecutiveDelimitersRule = "join";
            opts.LeadingDelimitersRule = "ignore";
            opts = setvaropts(opts, "Var10", "WhitespaceRule", "preserve");
            opts = setvaropts(opts, "Var10", "EmptyFieldRule", "auto");
            opts = setvaropts(opts, "X_PT", "TrimNonNumeric", true);
            opts = setvaropts(opts, "X_PT", "ThousandsSeparator", ",");
            temW = readtable(strcat(arps_dir,'weather_',num2str(t),'.txt'), opts);
            temW = table2array(temW);
            clear opts
        elseif(strcmp(b,'Volumes'))%my laptop
            opts = delimitedTextImportOptions("NumVariables", 10);
            opts.DataLines = [2, Inf];
            opts.Delimiter = " ";
            opts.VariableNames = ["X_PT", "Y_PT", "LATITUDE", "LONGITUDE", "TEMPA", "TEMPS", "RELH", "WSPD", "WDIR", "Var10"];
            opts.SelectedVariableNames = ["X_PT", "Y_PT", "LATITUDE", "LONGITUDE", "TEMPA", "TEMPS", "RELH", "WSPD", "WDIR"];
            opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "string"];
            opts.ExtraColumnsRule = "ignore";
            opts.EmptyLineRule = "read";
            opts.ConsecutiveDelimitersRule = "join";
            opts.LeadingDelimitersRule = "ignore";
            opts = setvaropts(opts, "Var10", "WhitespaceRule", "preserve");
            opts = setvaropts(opts, "Var10", "EmptyFieldRule", "auto");
            opts = setvaropts(opts, "X_PT", "TrimNonNumeric", true);
            opts = setvaropts(opts, "X_PT", "ThousandsSeparator", ",");
            temW = readtable(strcat(arps_dir,'weather_',num2str(t),'.txt'), opts);
            temW = table2array(temW);
            clear opts
        end
        for r=1:size(temW,1)
            if(~isnan(temW(r,2)))
                DEVS.WSPD(temW(r,2),temW(r,1),c)=temW(r,8);
                DEVS.WDIR(temW(r,2),temW(r,1),c)=temW(r,9);
            end
        end
        c=c+1;
    end
    [DEVS.U6,DEVS.V6]=sd2uv(DEVS.WSPD,DEVS.WDIR);
    DEVS = rmfield(DEVS,'WSPD');DEVS = rmfield(DEVS,'WDIR');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %6. Read in variables from ARPS history files: ARPS.HS, ARPS.U6, ARPS.V6
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    rst=10620;
    ren=10620+7200;
    rin=60;
    NSTYP=2;NZS=2;
    NT=(ren-rst)/rin+1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ARPS.PTSFLX=nan(NY,NX,NT);
    ARPS.U6=nan(NY,NX,NT);ARPS.V6=nan(NY,NX,NT);
    %Loop over files.
    c=1;
    for t=rst:rin:ren
        disp(num2str(t))
        if(t<10)
            numsec=strcat('00000',num2str(t));
        elseif(t<100)
            numsec=strcat('0000',num2str(t));
        elseif(t<1000)
            numsec=strcat('000',num2str(t));
        elseif(t<10000)
            numsec=strcat('00',num2str(t));
        elseif(t<100000)
            numsec=strcat('0',num2str(t));
        else
            numsec=strcat('',num2str(t));
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %mac server:
        if(c==1);temZP = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/zp','Index',{[1 1 1],[1 1 1],[NZ NY NX]});end
        temU = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/u','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
        temV = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/v','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
        temPTSFLX = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/ptsflx','Index',{[1 1],[1 1],[NY NX]});
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Unstagger the grid.
        %PT and QV (and thus, PTV) are defined at the *center* of the grid box
        %compute unstaggered U-wind.
        UWNDu=nan(NZ,NY,NX);
        for i=1:NX-1
            UWNDu(:,:,i)=0.5*(temU(:,:,i)+temU(:,:,i+1));
        end
        UWNDu(:,:,NX)=temU(:,:,NX);%Same as ARPS unstaggering (in gradsio3d.f90)
        %compute unstaggered V-wind.
        VWNDu=nan(NZ,NY,NX);
        for j=1:NY-1
            VWNDu(:,j,:)=0.5*(temV(:,j,:)+temV(:,j+1,:));
        end
        VWNDu(:,NY,:)=temV(:,NY,:);%Same as ARPS unstaggering (in gradsio3d.f90)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Linear vertical interpolation of U,V to 6 m AGL.
        ARPS.U6(:,:,c)=0.5*(UWNDu(4,:,:)+UWNDu(5,:,:));clear UWNDu temU
        ARPS.V6(:,:,c)=0.5*(VWNDu(4,:,:)+VWNDu(5,:,:));clear VWNDu temV
        %Convert sensible heat flux units to W/m2, and make positive
        %upward.  ARPS native units are K*kg*m-2*s-1.
        cp=1004.67;%Specific heat of air (constant pressure); Stull (1988).
        ARPS.HS(:,:,c)=temPTSFLX*(-cp);clear cp
        c=c+1;
    end
    clear temPTSFLX temW temZP
    %%%Mask values less than 500 W/m2
    ARPS.HS(ARPS.HS<500)=NaN;
    %%%Multiply wind components by factor of 0.7 to be consistent with what is
    %%%done in arps.f90:
    % !---------------------------------------------------------------------------
    % ! Begin MTK 7/9/24
    % !---------------------------------------------------------------------------
    % !Multiply wind components by factor of 0.7, to yield wind speeds more in line
    % !with what was observed at the control tower [see Fig. 5 in Kiefer et al.
    % !(2022)]
    % temfloat3 = 0.7 * temfloat3
    % temfloat4 = 0.7 * temfloat4
    % !---------------------------------------------------------------------------
    % ! End MTK 7/9/24
    % !---------------------------------------------------------------------------
    ARPS.U6 = 0.7 * ARPS.U6;
    ARPS.V6 = 0.7 * ARPS.V6;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Define plotting bounds.
    ss_s=39.9114;
    ss_n=39.9173;
    ss_w=-74.5980;
    ss_e=-74.5911;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Define SLEF_2019 burn block perimeter.  Coordinates obtained from
    %"SERDP_large_scale_finalconfig.kmz", created by reading in shapefiles in
    %"SERDP_large_scale_spatial_data" directory using Google Earth Pro.
    BBLOCK=[-74.59478356040911,39.91677961082065;
        -74.59457115429819,39.91666779013654;
        -74.59452728724368,39.91664469633786;
        -74.59286380514664,39.91576893296700;
        -74.59263123998637,39.91562093317582;
        -74.59245131082572,39.91549081265011;
        -74.59224486943265,39.91534151815286;
        -74.59219733516515,39.91532320076401;
        -74.59220090920040,39.91531563310795;
        -74.59226324345622,39.91515563632733;
        -74.59226872834930,39.91514155791109;
        -74.59214949636225,39.91441498988535;
        -74.59219582838328,39.91346853709193;
        -74.59214355304454,39.91285929948788;
        -74.59216447894477,39.91207669491202;
        -74.59229373771606,39.91177835865396;
        -74.59263215981019,39.91219631622555;
        -74.59274576883568,39.91224751218507;
        -74.59285219967526,39.91240270661773;
        -74.59304415190459,39.91256071757653;
        -74.59315115540666,39.91270613079016;
        -74.59355873977640,39.91290903402670;
        -74.59369321735265,39.91304491409569;
        -74.59384034732788,39.91319357752478;
        -74.59416305867801,39.91340705176342;
        -74.59478498058890,39.91398541412467;
        -74.59524827802721,39.91431262730358;
        -74.59556009467440,39.91463349654136;
        -74.59562972863725,39.91470515157850;
        -74.59570199754519,39.91477951783602;
        -74.59578102670557,39.91486084030277;
        -74.59595800804065,39.91499069234046;
        -74.59620860952352,39.91515513539763;
        -74.59653805978076,39.91547622531272;
        -74.59693773710107,39.91580114344045;
        -74.59701890627328,39.91604128965496;
        -74.59715001203868,39.91627435723552;
        -74.59673555383418,39.91634775951010;
        -74.59478574055204,39.91669305754058;
        -74.59478356040911,39.91677961082065];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Define time vector.
    tvec=linspace(datenum('2019-03-13 14:57'),datenum('2019-03-13 16:57'),121);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %SERDP FireTracker data (source: Eric Mueller)
    if(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server
        ft_dir='/dska/mtkiefer/SERDP/CASES/SLEF_2019/FIRE/';
    elseif(strcmp(b,'Volumes'))%my laptop
        ft_dir='/Volumes/SAMSUNG/SERDP_macbookair/CASES/SLEF_2019/FIRE/';
    end   
    load(strcat(ft_dir,'SLEF-0030_firetracker.mat'),'FT')
    %fire-tracker lat/lon
    FT.coords=[-74.59302562517473,39.91260969247371
        -74.59285623156497,39.91294520952201
        -74.59220477006532,39.91291866727881
        -74.59223814646226,39.91261978310636
        -74.59225408500302,39.91178624385721
        -74.5921904678668,39.91332438406685
        -74.59215272301118,39.9139710884039
        -74.5921678513252,39.91444916506787
        -74.59222480044716,39.91473421063718
        -74.59218991575436,39.91502862745747
        -74.59224791889537,39.9153352127505
        -74.59285880910143,39.91521927961341
        -74.5927624848003,39.91566017278522
        -74.59316403865088,39.91595396235233
        -74.59663344823109,39.91559864747305
        -74.59664317673342,39.91585178848273
        -74.59621448207342,39.91553260795672
        -74.59619261797937,39.91523690482378
        -74.59576380437905,39.9148964289751
        -74.59540495767475,39.91458047289041
        -74.59504852060834,39.91421996275075
        -74.5946734414629,39.91389433529815
        -74.59424899090126,39.91356995914325
        -74.5938578678564,39.91328627630555
        -74.59343055288149,39.91287478387966
        -74.59330813293407,39.91327197215425
        -74.59285307505104,39.91326948248405
        -74.59281397587094,39.91392129168687
        -74.59283001841922,39.91425260826344
        -74.59278247091429,39.91458234581287
        -74.59283189951108,39.91493784438832
        -74.59364825268571,39.91586342902342
        -74.59414301223109,39.91583684749192
        -74.59448543726246,39.91585344974929
        -74.59491153200001,39.91586394210832
        -74.59534453121194,39.91585512224363
        -74.59619629521087,39.91584645511166
        -74.59581455549815,39.91553547084972
        -74.59535063140528,39.91556416158762
        -74.59455416813945,39.91554468209529
        -74.59409086718287,39.91552589442343
        -74.59367617710112,39.91554341377928
        -74.59326446588138,39.91553567513164
        -74.59327025067113,39.91522282573548
        -74.59366947671846,39.91523222474945
        -74.59411057614041,39.91520032390732
        -74.59456426692054,39.91525115480713
        -74.59496790332406,39.9151948692189
        -74.59536110245604,39.91522252899478
        -74.595784051723,39.91520752349585
        -74.59538969803832,39.91490207885579
        -74.59495737626386,39.91489046457853
        -74.59456870739213,39.91488910631539
        -74.59410812827231,39.91487641000862
        -74.5936475493381,39.91486371187327
        -74.59324435439783,39.91488446837787
        -74.59323186095614,39.91456295793844
        -74.5936351195239,39.91453111664012
        -74.59409556535606,39.91456598446871
        -74.59455633847952,39.91454542633223
        -74.59495940047994,39.91454683490941
        -74.59454383895013,39.91422391598489
        -74.59409779016025,39.91418910033345
        -74.59366600581524,39.91417650306522
        -74.59324848548745,39.9141861243376
        -74.59323592655264,39.91387569868773
        -74.59368230035182,39.91385509350449
        -74.5940996223205,39.9138787251447];





    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Read in SERDP plant area density data (source: Kiefer et al.
    %2022, GMD)
    if(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server
        cn_dir='/dska/mtkiefer/SERDP/CASES/SLEF_2019/obs_data/LiDAR_CANOPY/PRODUCTION_CANOPY_PREP_vgrid_stretch_from_84m/';
    elseif(strcmp(b,'Volumes'))%my laptop
        cn_dir='/Volumes/SAMSUNG/SERDP_macbookair/CASES/SLEF_2019/obs_data/LiDAR_CANOPY/PRODUCTION_CANOPY_PREP_vgrid_stretch_from_84m/';
    end   
    tem=load(strcat(cn_dir,'SLEF-0030_arpscanopyfinal.mat'));
    CANOPY=tem.SLEF_30m;clear tem
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Save workspace for quick access.
    save(strcat(mat_dir,'/forpub/acfire_figures_forpub.mat'),'ARPS','DEVS',...
        'BBLOCK','FT','CANOPY','ss_n','ss_s','ss_w','ss_e','tvec')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Housekeeping
    clear arps_dir arps_head b c cmap devs_dir devsconfig_dir grid_dir i j 
    clear mat_dir NSTYP NT numsec NWRITEPTS_X NWRITEPTS_Y NX NY NZ NZS 
    clear plot_dir r ren RESRATX RESRATY rin rst t te ti ts ft_dir cn_dir
else
    load(strcat(mat_dir,'/forpub/acfire_figures_forpub.mat')) %REQUISITE MAT FILE FOUND
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PLOTTING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Introduce small adjustment to latitude and longitude arrays to correct for
%a ~5 m misalignment of the burn block in DEVS-FIRE.
adjY=5e-5;adjX=6.5e-5;
DEVS.LATS=DEVS.LATS-adjY;DEVS.LAT2=DEVS.LAT2-adjY;
DEVS.LONS=DEVS.LONS-adjX;DEVS.LON2=DEVS.LON2-adjX;
ARPS.LATS=ARPS.LATS-adjY;ARPS.LAT2=ARPS.LAT2-adjY;
ARPS.LONS=ARPS.LONS-adjX;ARPS.LON2=ARPS.LON2-adjX;

%%%FIGURE 1%%%
% Summary of PAD dataset implemented in ARPS-CANOPY.  
cmap=brewermap(16,'YlGn');%cmap(11,:)=[0.88 0.94 1];cmap(12,:)=[0.95 0.95 1];
close all
%%READ IN NECESSARY SHAPEFILES
%% Common variables
FontSize = 12;
FontName = 'MyriadPro-Regular'; % or choose any other font
doExportPlot = true;
%% start plot
figure_width = 6.5;
figure_height = 9;
% --- setup plot windows
figuresVisible = 'on'; % 'off' for non displayed plots (will still be exported)
hfig = figure(1); clf;
set(hfig,'Visible', figuresVisible)
set(hfig, 'units', 'inches', 'pos', [0 0 figure_width figure_height])
set(hfig, 'PaperPositionMode', 'auto');
set(hfig, 'Renderer','Painters');
set(hfig, 'Color', [1 1 1]); % Sets figure background
%set(gca, 'Color', [1 1 1]); % Sets axes background
h=usamap([ss_s+8e-5 ss_n-8e-5], [ss_w+1.75e-4 ss_e]);setm(h,'FontSize',10,'FFaceColor', [0.8 0.8 0.8],'LabelFormat','none','GLineWidth',1,'GColor','k');
hold on
%ARPS-CANOPY PLANT AREA DENSITY (PAI)
VAReps=squeeze(CANOPY.CLAIr_interp(1,:,:));VAReps(VAReps<1e-3)=NaN;
h1=pcolorm(double(CANOPY.LAT2),double(CANOPY.LON2),VAReps);h1.EdgeColor='none';
h1.FaceAlpha=1.0;
caxis([0 1.6]);colormap(cmap);
cb=colorbar('SouthOutside');
set(cb,'YTick',[0:0.5:1.5]);set(cb,'YTickLabel',[0:0.5:1.5]);
cb.FontSize=12;set(get(cb,'title'),'string','PAI','FontWeight','bold','FontAngle','normal','FontSize',12);
set(get(cb,'title'),'Position',[181.2500 8.25 0]);
set(cb,'Visible','on');
grid on
set(gca, 'layer', 'top','GridAlpha',0.35);
set(gca,'FontSize',12)
%ARPS-CANOPY CANOPY HEIGHT [METERS]
Zcan = [0:2:30];
VAReps=NaN*CANOPY.KCANOPEEr_interp;
for j=1:size(CANOPY.KCANOPEEr_interp,1)
    for i=1:size(CANOPY.KCANOPEEr_interp,2)
        if(~isnan(CANOPY.KCANOPEEr_interp(j,i)))
            VAReps(j,i)=Zcan(CANOPY.KCANOPEEr_interp(j,i));
        end
    end
end
VAReps(squeeze(CANOPY.CLAIr_interp(1,:,:))<1e-3)=NaN;%If PAI is very small, canopy height is undefined.
[C,h]=contourm(double(CANOPY.XLAT),double(CANOPY.XLON),double(VAReps),...
     'LevelList',[10:2:30],...
     'Fill','off','LineStyle','-','LineWidth',1,'LineColor','k','ShowText','on');
t = clabelm(C,h,'LabelSpacing',5000);
set(t,"BackgroundColor",[0.9 0.9 0.9 0.5]);%9/12/25: 0.5 alpha is a compromise between Sharon's suggestion that white background be completely removed and my opinion that no background color makes labels hard to read.
set(t,"Color",'k')
set(t,"FontWeight","bold")
set(t,"FontSize",10)
set(t,"Margin",1)
geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor','none','EdgeColor','k','LineWidth',3)
% Create rectangle *behind* north arrow (to make it easier to read)
gs=patchm([39.91152 39.91152 39.9121 39.9121 39.91150],[-74.5917 -74.59112 -74.59112 -74.5917 -74.5917],'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
northarrow('latitude', 39.9116, 'longitude', -74.5914,'scaleratio',0.075,'linewidth',2);
% Create rectangle *behind* scale ruler (to make it easier to read)
gs=patchm([39.91160 39.91160 39.91200 39.91200 39.91160],[-74.59755 -74.59550 -74.59550 -74.59755 -74.59755],'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
scaleruler on
setm(handlem('scaleruler1'), ...
    'XLoc',-216.974,'YLoc',4.793700e6,...
    'MajorTick',0:90:90,...
    'MajorTickLength',6,...
    'MinorTick',0:30:30,'Units','m',...
    'MinorTickLength',3,...
    'FontWeight','bold',...
    'LineWidth',2,'FontSize',14,'Color','k')
pos=cb.Position;cb.Position=[0.125 0.175 0.775 0.0125];
%Generate figure in PNG format (Painters: vector, not raster).
%TO DO: VECTOR IMAGE WITH TRANSPARENCY.
set(gcf, 'PaperPositionMode', 'Manual');
set(gcf, 'PaperUnits', 'Inches');
set(gcf,'Renderer','Painters')
set(gcf, 'PaperPosition',[1 1 6.5 9]);%Control dimensions of figure
h=gcf;
print(h,'-dpng',strcat(plot_dir,'/forpub/figure1.png'),'-r300');
%Quick cropping
A = imread(strcat(plot_dir,'/forpub/figure1.png'));
A1=imcrop(A, [198.5 451.5 1576 1864]);
imwrite(A1,strcat(plot_dir,'/forpub/figure1.png'));


%%%FIGURE 2%%%
% Burn progression map
cmap=brewermap(11,'Paired');
close all
%%READ IN NECESSARY SHAPEFILES
%% Common variables
FontSize = 12;
FontName = 'MyriadPro-Regular'; % or choose any other font
doExportPlot = true;
%% start plot
figure_width = 6.5;
figure_height = 9;
% --- setup plot windows
figuresVisible = 'on'; % 'off' for non displayed plots (will still be exported)
hfig = figure(1); clf;
set(hfig,'Visible', figuresVisible)
set(hfig, 'units', 'inches', 'pos', [0 0 figure_width figure_height])
set(hfig, 'PaperPositionMode', 'auto');
set(hfig, 'Renderer','Painters');
set(hfig, 'Color', [1 1 1]); % Sets figure background
%set(gca, 'Color', [1 1 1]); % Sets axes background
h=usamap([ss_s+8e-5 ss_n-8e-5], [ss_w+1.75e-4 ss_e]);setm(h,'FontSize',10,'FFaceColor', 'none','LabelFormat','none','GLineWidth',1,'GColor','k');
geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor',[0.5 0.5 0.5],'EdgeColor','k','LineWidth',3)
hold on
for c=1:size(FT.coords,1)
    s=scatterm(FT.coords(c,2),FT.coords(c,1),'filled','k');
    s.Children.MarkerFaceColor=[0.2 0.2 0.2];
    s.Children.MarkerEdgeColor='none';
    s.Parent.Clipping='off';
    s.Children.Clipping='off';
end
tem=DEVS.BPRG(1:499,1:499)-1;tem(tem>7200)=NaN;
c=contourm(DEVS.LATS(1:499,1:499),DEVS.LONS(1:499,1:499),tem,...
     'LevelList',[601:600:6601],...
     'Fill','off','LineStyle','-','LineWidth',1,'ShowText','off');%1:499: contourm cannot handle "NaNs" at N and E lateral boundaries.
caxis([601 6601]);colormap(cmap);
% Create rectangle *behind* legend (to make it easier to read)
annotation(hfig,'rectangle',...
    [0.79 0.29 0.08 0.275],...
    'LineStyle','-',...
    'EdgeColor','k',...
    'FaceColor',[0.5 0.5 0.5])
%Annotate: add legend with contour labels to east of burn block (time EDT).
annotation(hfig,'textbox',[0.785188034188037 0.295 0.112247863247861 0.0246913580246913],'String','15:07','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(1,:));
annotation(hfig,'textbox',[0.785188034188037 0.320 0.112247863247861 0.0246913580246913],'String','15:17','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(2,:));
annotation(hfig,'textbox',[0.785188034188037 0.3450 0.112247863247861 0.0246913580246913],'String','15:27','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(3,:));
annotation(hfig,'textbox',[0.785188034188037 0.3700 0.112247863247861 0.0246913580246913],'String','15:37','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(4,:));
annotation(hfig,'textbox',[0.785188034188037 0.3950 0.112247863247861 0.0246913580246913],'String','15:47','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(5,:));
annotation(hfig,'textbox',[0.785188034188037 0.4200 0.112247863247861 0.0246913580246913],'String','15:57','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(6,:));
annotation(hfig,'textbox',[0.785188034188037 0.4450 0.112247863247861 0.0246913580246913],'String','16:07','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(7,:));
annotation(hfig,'textbox',[0.785188034188037 0.4700 0.112247863247861 0.0246913580246913],'String','16:17','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(8,:));
annotation(hfig,'textbox',[0.785188034188037 0.4950 0.112247863247861 0.0246913580246913],'String','16:27','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(9,:));
annotation(hfig,'textbox',[0.785188034188037 0.5200 0.112247863247861 0.0246913580246913],'String','16:37','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(10,:));
annotation(hfig,'textbox',[0.785188034188037 0.5450 0.112247863247861 0.0246913580246913],'String','16:47','FontSize',12,'FontWeight','bold','LineStyle','none','FitBoxToText','off','Color',cmap(11,:));
%
grid on
set(gca, 'layer', 'top','GridAlpha',0.35);
% Create rectangle *behind* north arrow (to make it easier to read)
%gs=patchm([39.91152 39.91152 39.9121 39.9121 39.91150],[-74.5917 -74.59112 -74.59112 -74.5917 -74.5917],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
northarrow('latitude', 39.9116, 'longitude', -74.5914,'scaleratio',0.075,'linewidth',2);
% Create rectangle *behind* scale ruler (to make it easier to read)
%gs=patchm([39.91160 39.91160 39.91200 39.91200 39.91160],[-74.59755 -74.59550 -74.59550 -74.59755 -74.59755],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
scaleruler on
setm(handlem('scaleruler1'), ...
    'XLoc',-216.974,'YLoc',4.793700e6,...
    'MajorTick',0:90:90,...
    'MajorTickLength',6,...
    'MinorTick',0:30:30,'Units','m',...
    'MinorTickLength',3,...
    'FontWeight','bold',...
    'LineWidth',2,'FontSize',14,'Color','k')
set(gca,'FontSize',12)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Firetracker overlay%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ts=(3600*(14+4))+(57*60)+0;%seconds since 00 UTC: 14:57:00 EDT.
%te=(3600*(16+4))+(57*60)+0;%seconds since 00 UTC: 16:57:00 EDT.
ti=600;%contour fire front position every 10 minutes (midpoints of colorbar ranges)
tem1=FT.timeCt;tem=datenum('2019-03-13 00:00','yyyy-mm-dd HH:MM')+(tem1/86400);clear tem1%Convert to datenum object
ts=datenum('2019-03-13 19:07','yyyy-mm-dd HH:MM');%15:07 EDT
te=datenum('2019-03-13 20:47','yyyy-mm-dd HH:MM');%16:47 EDT
c=1;
for levs=ts:(ti/86400):te%fractional day
    contourm(FT.latd,FT.lond,tem,levs,'LineWidth',1,'LineColor',cmap(c,:),'LineStyle','--')
    c=c+1;
end
%Depict outline of inset panels.
geoshow([39.9151 39.9160 39.9160 39.9151 39.9151],[-74.5953 -74.5953 -74.5946 -74.5946 -74.5953],'DisplayType','polygon','FaceColor','none','EdgeColor','k','LineWidth',2,'LineStyle',':')
geoshow([39.9139 39.9148 39.9148 39.9139 39.9139],[-74.5928 -74.5928 -74.5921 -74.5921 -74.5928],'DisplayType','polygon','FaceColor','none','EdgeColor','k','LineWidth',2,'LineStyle',':')
%//////Begin inset panel 1//////
axes('Position',[.175 .395 .18 .18])
box on
h=usamap([39.9151 39.9160], [-74.5953  -74.5946]);setm(h,'FontSize',10,'FFaceColor', 'none','LabelFormat','none','GLineWidth',1,'GColor','k');
geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor',[0.5 0.5 0.5],'EdgeColor','k','LineWidth',1)
hold on
for c=1:size(FT.coords,1)
    s=scatterm(FT.coords(c,2),FT.coords(c,1),'filled','k');
    s.Children.MarkerFaceColor=[0.2 0.2 0.2];
    s.Children.MarkerEdgeColor='none';
    s.Parent.Clipping='off';
    s.Children.Clipping='off';
end
tem=DEVS.BPRG(1:499,1:499)-1;tem(tem>7200)=NaN;
c=contourm(DEVS.LATS(1:499,1:499),DEVS.LONS(1:499,1:499),tem,...
     'LevelList',[601:600:6601],...
     'Fill','off','LineStyle','-','LineWidth',1,'ShowText','off');%1:499: contourm cannot handle "NaNs" at N and E lateral boundaries.
caxis([601 6601]);colormap(cmap);
set(gca,'FontSize',12)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Firetracker overlay%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ts=(3600*(14+4))+(57*60)+0;%seconds since 00 UTC: 14:57:00 EDT.
%te=(3600*(16+4))+(57*60)+0;%seconds since 00 UTC: 16:57:00 EDT.
ti=600;%contour fire front position every 10 minutes (midpoints of colorbar ranges)
tem1=FT.timeCt;tem=datenum('2019-03-13 00:00','yyyy-mm-dd HH:MM')+(tem1/86400);clear tem1%Convert to datenum object
ts=datenum('2019-03-13 19:07','yyyy-mm-dd HH:MM');%15:07 EDT
te=datenum('2019-03-13 20:47','yyyy-mm-dd HH:MM');%16:47 EDT
c=1;
for levs=ts:(ti/86400):te%fractional day
    contourm(FT.latd,FT.lond,tem,levs,'LineWidth',1,'LineColor',cmap(c,:),'LineStyle','--')
    c=c+1;
end
scaleruler on
setm(handlem('scaleruler1'), ...
    'XLoc',-10.9966,'YLoc',4.794115e6,...
    'MajorTick',0:30:30,...
    'MajorTickLength',6,...
    'MinorTick',0:10:10,'Units','m',...
    'MinorTickLength',3,...
    'FontWeight','bold',...
    'LineWidth',2,'FontSize',9,'Color','k')
set(gca,'FontSize',12)
uistack(h,'top')
%\\\\\End inset panel 1\\\\\
%//////Begin inset panel 2//////
axes('Position',[.3925 .225 .18 .18])
box on
h=usamap([39.9139 39.9148], [-74.5928 -74.5921]);setm(h,'FontSize',10,'FFaceColor', 'none','LabelFormat','none','GLineWidth',1,'GColor','k');
geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor',[0.5 0.5 0.5],'EdgeColor','k','LineWidth',1)
hold on
for c=1:size(FT.coords,1)
    s=scatterm(FT.coords(c,2),FT.coords(c,1),'filled','k');
    s.Children.MarkerFaceColor=[0.2 0.2 0.2];
    s.Children.MarkerEdgeColor='none';
    s.Parent.Clipping='off';
    s.Children.Clipping='off';
end
tem=DEVS.BPRG(1:499,1:499)-1;tem(tem>7200)=NaN;
c=contourm(DEVS.LATS(1:499,1:499),DEVS.LONS(1:499,1:499),tem,...
     'LevelList',[601:600:6601],...
     'Fill','off','LineStyle','-','LineWidth',1,'ShowText','off');%1:499: contourm cannot handle "NaNs" at N and E lateral boundaries.
caxis([601 6601]);colormap(cmap);
set(gca,'FontSize',12)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Firetracker overlay%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ts=(3600*(14+4))+(57*60)+0;%seconds since 00 UTC: 14:57:00 EDT.
%te=(3600*(16+4))+(57*60)+0;%seconds since 00 UTC: 16:57:00 EDT.
ti=600;%contour fire front position every 10 minutes (midpoints of colorbar ranges)
tem1=FT.timeCt;tem=datenum('2019-03-13 00:00','yyyy-mm-dd HH:MM')+(tem1/86400);clear tem1%Convert to datenum object
ts=datenum('2019-03-13 19:07','yyyy-mm-dd HH:MM');%15:07 EDT
te=datenum('2019-03-13 20:47','yyyy-mm-dd HH:MM');%16:47 EDT
c=1;
for levs=ts:(ti/86400):te%fractional day
    contourm(FT.latd,FT.lond,tem,levs,'LineWidth',1,'LineColor',cmap(c,:),'LineStyle','--')
    c=c+1;
end
scaleruler on
setm(handlem('scaleruler1'), ...
    'XLoc',-12.4432,'YLoc',4.793950e6,...
    'MajorTick',0:30:30,...
    'MajorTickLength',6,...
    'MinorTick',0:10:10,'Units','m',...
    'MinorTickLength',3,...
    'FontWeight','bold',...
    'LineWidth',2,'FontSize',9,'Color','k')
set(gca,'FontSize',12)
uistack(h,'top')
%\\\\\End inset panel 2\\\\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Generate figure in PNG format (Painters: vector, not raster).
%TO DO: VECTOR IMAGE WITH TRANSPARENCY.
set(gcf, 'PaperPositionMode', 'Manual');
set(gcf, 'PaperUnits', 'Inches');
set(gcf,'Renderer','Painters')
set(gcf, 'PaperPosition',[1 1 6.5 9]);%Control dimensions of figure
h=gcf;
print(h,'-dpng',strcat(plot_dir,'/forpub/figure2.png'),'-r300');
%Quick cropping
A = imread(strcat(plot_dir,'/forpub/figure2.png'));
A1=imcrop(A, [251.5 455.5 1520 1696]);
imwrite(A1,strcat(plot_dir,'/forpub/figure2.png'));


%%%FIGURES 3-8%%%
% Heat flux / wind field overlay plots
%STRATEGY: PLOT HEAT FLUX - WIND FIELD OVERLAYS EVERY TWENTY MINUTES:
%FIG 3: 1507 EDT (t=11)
%FIG 4: 1527 EDT (t=31)
%FIG 5: 1547 EDT (t=51)
%FIG 6: 1607 EDT (t=71)
%FIG 7: 1627 EDT (t=91)
%FIG 8: 1647 EDT (t=111)
c=3;
for t=[11:20:111]
    close all
    %%READ IN NECESSARY SHAPEFILES
    %% Common variables
    FontSize = 12;
    FontName = 'MyriadPro-Regular'; % or choose any other font
    doExportPlot = true;
    %% start plot
    figure_width = 6.5;
    figure_height = 9;
    % --- setup plot windows
    figuresVisible = 'on'; % 'off' for non displayed plots (will still be exported)
    hfig = figure(1); clf;
    set(hfig,'Visible', figuresVisible)
    set(hfig, 'units', 'inches', 'pos', [0 0 figure_width figure_height])
    set(hfig, 'PaperPositionMode', 'auto');
    set(hfig, 'Renderer','Painters');
    set(hfig, 'Color', [1 1 1]); % Sets figure background
    %set(gca, 'Color', [1 1 1]); % Sets axes background
    h=usamap([ss_s+8e-5 ss_n-8e-5], [ss_w+1.75e-4 ss_e]);setm(h,'FontSize',10,'FFaceColor', 'none','LabelFormat','none','GLineWidth',1,'GColor','k');
    hold on
    %DEVS-FIRE sensible heat flux
    h1=pcolorm(double(DEVS.LAT2),double(DEVS.LON2),squeeze(DEVS.HS(:,:,t))/1000.);h1.EdgeColor='none';%W/m2 --> kW/m2
    h1.FaceAlpha=1.0;
    %ARPS-CANOPY sensible heat flux
    h2=pcolorm(double(ARPS.LAT2),double(ARPS.LON2),squeeze(ARPS.HS(:,:,t))/1000.);h2.EdgeColor='none';%W/m2 --> kW/m2
    h2.FaceAlpha=0.5;
    caxis([0 165]);colormap(brewermap(33,'*Spectral'));
    cb=colorbar('SouthOutside');set(cb,'YTick',[0:10:160]);
    set(cb,'YTickLabel',[0.5,10:10:160]);set(cb,'TickLength',0);cb.Ruler.TickLabelRotation=0;
    cb.FontSize=12;set(get(cb,'title'),'string','Hs [kW m^{-2}]','FontWeight','bold','FontAngle','normal','FontSize',12);
    set(get(cb,'title'),'Position',[181.2500 8.25 0]);
    set(cb,'Visible','on');
    grid on
    tt=text(0.5,1.02,{[datestr(tvec(t),'dd-mmm-yyyy HH:MM:SS'),' EDT']},'Units','normalized','HorizontalAlignment','center');tt.FontSize=14;tt.FontWeight='bold';tt.Color='k';
    set(gca, 'layer', 'top','GridAlpha',0.35);
    set(gca,'FontSize',12)
    geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor','none','EdgeColor','k','LineWidth',2)
    %DEVS-FIRE wind field (plot every third grid point)
    q1=quivermc(double(DEVS.LATS),double(DEVS.LONS),squeeze(DEVS.U6(:,:,t)),squeeze(DEVS.V6(:,:,t)),'Color',[0.1 0.1 0.1],'LineWidth',0.5,'Density',33,'reference',3,'refvec','true','arrowstyle','tail');
    %ARPS-CANOPY wind field (plot every grid point)
    q2=quivermc(double(ARPS.LATS),double(ARPS.LONS),squeeze(ARPS.U6(:,:,t)),squeeze(ARPS.V6(:,:,t)),'Color',[0.0 0.1 0.9],'LineWidth',1.0,'Density',100,'reference',3,'refvec','true','arrowstyle','tail');
    % Create rectangle *behind* north arrow (to make it easier to read)
    gs=patchm([39.91152 39.91152 39.9121 39.9121 39.91150],[-74.5917 -74.59112 -74.59112 -74.5917 -74.5917],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
    n=northarrow('latitude', 39.9116, 'longitude', -74.5914,'scaleratio',0.075,'linewidth',2);n(1).ZData=2+0*n(1).XData;n(2).ZData=2+0*n(2).XData;%trick: move north arrow to top.
    % Create rectangle *behind* scale ruler (to make it easier to read)
    gs=patchm([39.91160 39.91160 39.91200 39.91200 39.91160],[-74.59755 -74.59550 -74.59550 -74.59755 -74.59755],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
    s=scaleruler;
    setm(s,'XLoc',-216.974,'YLoc',4.793700e6,'MajorTick',0:90:90,'MajorTickLength',6,'MinorTick',0:30:30,'Units','m','MinorTickLength',3,'FontWeight','bold','LineWidth',2,'FontSize',14,'Color','k');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Brute force approach to move scaleruler above patch (and, in turn, above quiver).
    for ss=[1 4 5 6];pos=s.Children(ss).Position;s.Children(ss).Position=[pos(1) pos(2) 2.5];clear pos;end
    for ss=[2 3 7 8];s.Children(ss).ZData=2.5+0*s.Children(ss).XData;end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cb.Position=[0.125 0.175 0.775 0.0125];
    %Generate figure in PNG format (Painters: vector, not raster).
    %TO DO: VECTOR IMAGE WITH TRANSPARENCY.
    set(gcf, 'PaperPositionMode', 'Manual');
    set(gcf, 'PaperUnits', 'Inches');
    set(gcf,'Renderer','Painters')
    set(gcf, 'PaperPosition',[1 1 6.5 9]);%Control dimensions of figure
    h=gcf;
    print(h,'-dpng',strcat(plot_dir,'/forpub/figure',num2str(c),'.png'),'-r300');
    %Quick cropping
    A = imread(strcat(plot_dir,'/forpub/figure',num2str(c),'.png'));
    A1=imcrop(A, [191.5 379.5 1596 1940]);
    imwrite(A1,strcat(plot_dir,'/forpub/figure',num2str(c),'.png'));
    c=c+1;
end




%%%SUPPLEMENTAL MATERIAL%%%
%MOVIE S1: EVERY MINUTE FROM 14:57 TO 16:57
for t=1:size(DEVS.HS,3)
    close all
    %%READ IN NECESSARY SHAPEFILES
    %% Common variables
    FontSize = 12;
    FontName = 'MyriadPro-Regular'; % or choose any other font
    doExportPlot = true;
    %% start plot
    figure_width = 6.5;
    figure_height = 9;
    % --- setup plot windows
    figuresVisible = 'on'; % 'off' for non displayed plots (will still be exported)
    hfig = figure(1); clf;
    set(hfig,'Visible', figuresVisible)
    set(hfig, 'units', 'inches', 'pos', [0 0 figure_width figure_height])
    set(hfig, 'PaperPositionMode', 'auto');
    set(hfig, 'Renderer','Painters');
    set(hfig, 'Color', [1 1 1]); % Sets figure background
    %set(gca, 'Color', [1 1 1]); % Sets axes background
    h=usamap([ss_s+8e-5 ss_n-8e-5], [ss_w+1.75e-4 ss_e]);setm(h,'FontSize',10,'FFaceColor', 'none','LabelFormat','none','GLineWidth',1,'GColor','k');
    hold on
    %DEVS-FIRE sensible heat flux
    h1=pcolorm(double(DEVS.LAT2),double(DEVS.LON2),squeeze(DEVS.HS(:,:,t))/1000.);h1.EdgeColor='none';%W/m2 --> kW/m2
    h1.FaceAlpha=1.0;
    %ARPS-CANOPY sensible heat flux
    h2=pcolorm(double(ARPS.LAT2),double(ARPS.LON2),squeeze(ARPS.HS(:,:,t))/1000.);h2.EdgeColor='none';%W/m2 --> kW/m2
    h2.FaceAlpha=0.5;
    caxis([0 165]);colormap(brewermap(33,'*Spectral'));
    cb=colorbar('SouthOutside');set(cb,'YTick',[0:10:160]);
    set(cb,'YTickLabel',[0.5,10:10:160]);set(cb,'TickLength',0);cb.Ruler.TickLabelRotation=0;
    cb.FontSize=12;set(get(cb,'title'),'string','Hs [kW m^{-2}]','FontWeight','bold','FontAngle','normal','FontSize',12);
    set(get(cb,'title'),'Position',[181.2500 8.25 0]);
    set(cb,'Visible','on');
    grid on
    tt=text(0.5,1.02,{[datestr(tvec(t),'dd-mmm-yyyy HH:MM:SS'),' EDT']},'Units','normalized','HorizontalAlignment','center');tt.FontSize=14;tt.FontWeight='bold';tt.Color='k';
    set(gca, 'layer', 'top','GridAlpha',0.35);
    set(gca,'FontSize',12)
    geoshow(BBLOCK(:,2),BBLOCK(:,1),'DisplayType','polygon','FaceColor','none','EdgeColor','k','LineWidth',2)
    %DEVS-FIRE wind field (plot every third grid point)
    q1=quivermc(double(DEVS.LATS),double(DEVS.LONS),squeeze(DEVS.U6(:,:,t)),squeeze(DEVS.V6(:,:,t)),'Color',[0.1 0.1 0.1],'LineWidth',0.5,'Density',33,'reference',3,'refvec','true','arrowstyle','tail');
    %ARPS-CANOPY wind field (plot every grid point)
    q2=quivermc(double(ARPS.LATS),double(ARPS.LONS),squeeze(ARPS.U6(:,:,t)),squeeze(ARPS.V6(:,:,t)),'Color',[0.0 0.1 0.9],'LineWidth',1.0,'Density',100,'reference',3,'refvec','true','arrowstyle','tail');
    % Create rectangle *behind* north arrow (to make it easier to read)
    gs=patchm([39.91152 39.91152 39.9121 39.9121 39.91150],[-74.5917 -74.59112 -74.59112 -74.5917 -74.5917],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
    n=northarrow('latitude', 39.9116, 'longitude', -74.5914,'scaleratio',0.075,'linewidth',2);n(1).ZData=2+0*n(1).XData;n(2).ZData=2+0*n(2).XData;%trick: move north arrow to top.
    % Create rectangle *behind* scale ruler (to make it easier to read)
    gs=patchm([39.91160 39.91160 39.91200 39.91200 39.91160],[-74.59755 -74.59550 -74.59550 -74.59755 -74.59755],1,'w');gs.EdgeColor='none';gs.FaceAlpha=0.75;
    s=scaleruler;
    setm(s,'XLoc',-216.974,'YLoc',4.793700e6,'MajorTick',0:90:90,'MajorTickLength',6,'MinorTick',0:30:30,'Units','m','MinorTickLength',3,'FontWeight','bold','LineWidth',2,'FontSize',14,'Color','k');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Brute force approach to move scaleruler above patch (and, in turn, above quiver).
    for ss=[1 4 5 6];pos=s.Children(ss).Position;s.Children(ss).Position=[pos(1) pos(2) 2.5];clear pos;end
    for ss=[2 3 7 8];s.Children(ss).ZData=2.5+0*s.Children(ss).XData;end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cb.Position=[0.125 0.175 0.775 0.0125];
    %Generate figure in PNG format (Painters: vector, not raster).
    %TO DO: VECTOR IMAGE WITH TRANSPARENCY.
    set(gcf, 'PaperPositionMode', 'Manual');
    set(gcf, 'PaperUnits', 'Inches');
    set(gcf,'Renderer','Painters')
    set(gcf, 'PaperPosition',[1 1 6.5 9]);%Control dimensions of figure
    h=gcf;
    print(h,'-dpng',strcat(plot_dir,'/forpub/movieS1/movie_S1_t',sprintf('%03d',t),'.png'),'-r300');
    %Quick cropping
    A = imread(strcat(plot_dir,'/forpub/movieS1/movie_S1_t',sprintf('%03d',t),'.png'));
    A1=imcrop(A, [191.5 379.5 1596 1940]);
    imwrite(A1,strcat(plot_dir,'/forpub/movieS1/movie_S1_t',sprintf('%03d',t),'.png'));
end
